Only single-component vaccines should be used for the birth dose and doses administered before age 6 weeks. (Source: CDC)
Combined hepatitis A and B vaccine is licensed for children from 1 year old. (Source: VNVC)
4 in 1, 5 in 1, and 6 in 1 vaccines are licensed for children aged 6 weeks through 6 years. (Source: CDC)
According to the Vietnamese vaccination schedule, 4 in 1, 5 in 1, and 6 in 1 vaccines should be administered in months 2, 3, 4. (Source: Tiem Chung Mo Rong)
library(tidyverse)
library(readxl)
library(writexl)
library(janitor)
library(Hmisc)
library(arrow)
library(gtsummary)
library(plotly)
library(lubridate)
library(data.table)
library(summarytools)# Import and review data
data <- read_parquet(here::here("data/hepb_data_long.parquet")) %>% collect()## Warning: Potentially unsafe or invalid elements have been discarded from R metadata.
## ℹ Type: "externalptr"
## → If you trust the source, you can set `options(arrow.unsafe_metadata = TRUE)` to preserve them.print(data)## pid vacname vacdate vactype
## <char> <char> <Date> <char>
## 1: 221194720170013 Quinvaxem 2017-07-20 TCCD
## 2: 221031520170032 Quinvaxem 2017-07-20 TCCD
## 3: 221031520170019 Quinvaxem 2017-07-20 TCCD
## 4: 221031520170025 Quinvaxem 2017-07-20 TCCD
## 5: 221131920170027 Quinvaxem 2017-07-19 TCCD
## ---
## 40421967: 217182920220004 SII 2022-05-17 TCMR
## 40421968: 217202720220001 SII 2022-05-18 TCMR
## 40421969: 217011420220019 Hep B vaccine for newborn 2022-05-19 TCMR
## 40421970: 217200520220010 SII 2022-05-17 TCMR
## 40421971: 217151520210070 SII 2022-05-16 TCMR
## vacplace vacplace_type province_reg
## <char> <char> <char>
## 1: TYT Đồng Việt TCMR Bắc Giang
## 2: TYT Tam Tiến TCMR Bắc Giang
## 3: TYT Tam Tiến TCMR Bắc Giang
## 4: TYT Tam Tiến TCMR Bắc Giang
## 5: TYT Vĩnh Khương-đã sát nhập TCMR Bắc Giang
## ---
## 40421967: TYT Tu Vũ TCMR Phú Thọ
## 40421968: TYT Xuân Sơn TCMR Phú Thọ
## 40421969: Bệnh viện Sản Nhi tỉnh Phú Thọ BV Phú Thọ
## 40421970: TYT Thu Cúc TCMR Phú Thọ
## 40421971: TYT Xuân Thủy TCMR Phú Thọ
## district_reg commune_reg sex dob ethnic
## <char> <char> <char> <Date> <char>
## 1: Yên Dũng Đồng Việt nam 2017-02-03 Kinh
## 2: Yên Thế Tam Tiến nu 2017-04-21 Kinh
## 3: Yên Thế Tam Tiến nam 2017-03-09 Kinh
## 4: Yên Thế Tam Tiến nam 2017-03-30 Kinh
## 5: Sơn Động Vĩnh An nam 2017-04-14 Kinh
## ---
## 40421967: Thanh Thủy Tu Vũ nam 2022-01-12 Kinh
## 40421968: Tân Sơn Xuân Sơn nu 2022-01-13 Mường
## 40421969: Thành phố Việt Trì Bến Gót nam 2022-05-19 Kinh
## 40421970: Tân Sơn Thu Cúc nam 2022-01-20 Mường
## 40421971: Yên Lập Xuân Thủy nu 2021-11-24 Kinhstr(data)## Classes 'data.table' and 'data.frame': 40421971 obs. of 12 variables:
## $ pid : chr "221194720170013" "221031520170032" "221031520170019" "221031520170025" ...
## $ vacname : chr "Quinvaxem" "Quinvaxem" "Quinvaxem" "Quinvaxem" ...
## $ vacdate : Date, format: "2017-07-20" "2017-07-20" ...
## $ vactype : chr "TCCD" "TCCD" "TCCD" "TCCD" ...
## $ vacplace : chr "TYT Đồng Việt" "TYT Tam Tiến" "TYT Tam Tiến" "TYT Tam Tiến" ...
## $ vacplace_type: chr "TCMR" "TCMR" "TCMR" "TCMR" ...
## $ province_reg : chr "Bắc Giang" "Bắc Giang" "Bắc Giang" "Bắc Giang" ...
## $ district_reg : chr "Yên Dũng" "Yên Thế" "Yên Thế" "Yên Thế" ...
## $ commune_reg : chr "Đồng Việt" "Tam Tiến" "Tam Tiến" "Tam Tiến" ...
## $ sex : chr "nam" "nu" "nam" "nam" ...
## $ dob : Date, format: "2017-02-03" "2017-04-21" ...
## $ ethnic : chr "Kinh" "Kinh" "Kinh" "Kinh" ...# Remove Boostrix, not HepB vaccine
data <- data %>% filter(vacname != "Boostrix")
# Check the number of entries not within the data collection time interval
data %>% filter(year(vacdate) < 2014 | year(vacdate) > 2022) %>% nrow()## [1] 622# Select only injection correct entries from 2014 to 2022
data <- data %>% filter(year(vacdate) >= 2014 & year(vacdate) <= 2022)#check missing values
missing_values <- is.na(data)
missing_counts <- colSums(is.na(data))
rows_with_missing <- data[!complete.cases(data), ]
print(rows_with_missing)## pid vacname vacdate vactype
## <char> <char> <Date> <char>
## 1: 221192720170048 Quinvaxem 2017-07-20 TCCD
## 2: 221011820190065 Hexaxim 2019-06-23 TCMR
## 3: 221153120190139 Hep B vaccine for newborn 2019-06-05 TCCD
## 4: 103051120180363 Hexaxim 2019-06-27 TCMR
## 5: 805011920210311 Hep B vaccine for newborn 2021-07-05 TCMR
## ---
## 1735181: 217110120220004 Infanrix Hexa 2022-03-19 TCMR
## 1735182: 217134320210195 SII 2022-03-16 TCCD
## 1735183: 217113720210100 Hexaxim 2021-12-15 TCDV
## 1735184: 217181720210236 Infanrix Hexa 2022-02-15 TCMR
## 1735185: 217010720210145 Hexaxim 2021-11-30 TCMR
## vacplace vacplace_type province_reg
## <char> <char> <char>
## 1: <NA> <NA> Bắc Giang
## 2: Trung tâm tiêm chủng VNVC Skyline Văn Quán DV Bắc Giang
## 3: <NA> <NA> Bắc Giang
## 4: Trung tâm tiêm chủng VNVC 180 Trường Chinh DV Bắc Giang
## 5: <NA> <NA> An Giang
## ---
## 1735181: <NA> <NA> Phú Thọ
## 1735182: <NA> <NA> Phú Thọ
## 1735183: <NA> Dia diem khac Phú Thọ
## 1735184: <NA> <NA> Phú Thọ
## 1735185: <NA> <NA> Phú Thọ
## district_reg commune_reg sex dob ethnic
## <char> <char> <char> <Date> <char>
## 1: Yên Dũng Tiền Phong nu 2017-03-18 Kinh
## 2: <NA> Dĩnh Trì nam 2019-01-07 Kinh
## 3: Lục Nam Nghĩa Phương nam 2019-06-04 Kinh
## 4: <NA> Trần Nguyên Hãn nu 2018-11-03 Kinh
## 5: Thành phố Long Xuyên Mỹ Hòa Hưng nu 2021-07-05 Kinh
## ---
## 1735181: Phù Ninh Thị trấn Phong Châu nu 2022-01-04 Kinh
## 1735182: Cẩm Khê Tạ Xá nam 2021-12-12 Kinh
## 1735183: Phù Ninh Bình Phú nam 2021-08-23 Kinh
## 1735184: Thanh Thủy Sơn Thủy nam 2021-10-05 Kinh
## 1735185: Thành phố Việt Trì Gia Cẩm nu 2021-07-04 Kinhyear_dob and year_vac are the years in which the child was born and vaccinated, respectively.
# Mutate columns for year of birth and year of vaccination
data <- data %>% mutate(year_dob = year(dob), year_vac = year(vacdate))Representing the age in days at which the child received the vaccine.
data <- data %>% mutate(vac_age = difftime(vacdate, dob, units = "days"))Representing the month after birth that the child receives the vaccine.
data <- data %>% mutate(vac_month = case_when(
vac_age < 0 ~ NA,
vac_age <= 1 ~ -2,
vac_age <= 7 & vac_age > 1 ~ -1,
vac_age < 30 & vac_age > 7 ~ 0,
vac_age >= 30 & vac_age < 60 ~ 1,
vac_age >= 60 & vac_age < 90 ~ 2,
vac_age >= 90 & vac_age < 120 ~ 3,
vac_age >= 120 & vac_age < 150 ~ 4,
vac_age >= 150 & vac_age < 180 ~ 5,
vac_age >= 180 & vac_age < 210 ~ 6,
vac_age >= 210 & vac_age < 240 ~ 7,
vac_age >= 240 & vac_age < 270 ~ 8,
vac_age >= 270 & vac_age < 300 ~ 9,
vac_age >= 300 & vac_age < 330 ~ 10,
vac_age >= 330 ~ 11
))Vietnam is officially divided into six regions (source1,source2, source3.)
data <- data %>%
mutate(region = case_when(
province_reg %in% c("Hà Nội", "Vĩnh Phúc", "Bắc Ninh", "Quảng Ninh", "Hải Dương", "Hải Phòng", "Hưng Yên", "Thái Bình", "Hà Nam", "Nam Định", "Ninh Bình") ~ "RRD",
province_reg %in% c("Hà Giang", "Cao Bằng", "Bắc Kạn", "Tuyên Quang", "Lào Cai", "Yên Bái", "Thái Nguyên", "Lạng Sơn", "Bắc Giang", "Phú Thọ", "Điện Biên", "Lai Châu", "Sơn La", "Hòa Bình") ~ "NMM",
province_reg %in% c("Thanh Hóa", "Nghệ An", "Hà Tĩnh", "Quảng Bình", "Quảng Trị", "Thừa Thiên Huế", "Đà Nẵng", "Quảng Nam", "Quảng Ngãi", "Bình Định", "Phú Yên", "Khánh Hòa", "Ninh Thuận", "Bình Thuận") ~ "NCC",
province_reg %in% c("Kon Tum", "Gia Lai", "Đắk Lắk", "Đắk Nông", "Lâm Đồng") ~ "CHL",
province_reg %in% c("Bình Phước", "Tây Ninh", "Bình Dương", "Đồng Nai", "Bà Rịa - Vũng Tàu", "Thành phố Hồ Chí Minh") ~ "SE",
province_reg %in% c("Long An", "Tiền Giang", "Bến Tre", "Trà Vinh", "Vĩnh Long", "Đồng Tháp", "An Giang", "Kiên Giang", "Cần Thơ", "Hậu Giang", "Sóc Trăng", "Bạc Liêu", "Cà Mau") ~ "MKD",
.default = NULL
)) %>%
mutate(region = factor(region, levels = c("RRD", "NMM", "NCC", "CHL", "SE", "MKD")))# Classify vacname and create a new variable vacgroup
data <- data %>%
mutate(vacgroup = case_when(
vacname %in% c("ENGERIX-B", "Euvax B", "Gene-Hbvax", "H-B-VAX II", "HBVaxPRO", "r-HBvax", "Heberbiovac HB", "Hep B vaccine", "Hep B vaccine for newborn", "SCI-B-VAC", "Hepavax-Gene TF") ~ "single",
vacname == "Twinrix" ~ "hepab",
vacname == "TRITANRIX-HB" ~ "4_in_1",
vacname %in% c("Quinvaxem", "SII") ~ "5_in_1",
vacname %in% c("Hexaxim", "Hexavac", "Infanrix Hexa") ~ "6_in_1",
.default = NULL # default case
))data <- data %>% mutate(vacname2 = case_when(
vacname %in% c("ENGERIX-B", "Euvax B", "Gene-Hbvax", "Heberbiovac HB", "Hep B vaccine", "Hepavax-Gene TF") & vac_age < 30 ~ "possible_hep_b_vaccine_for_newborn",
.default = vacname))
data %>% filter(vacname2 == "possible_hep_b_vaccine_for_newborn") %>% nrow()## [1] 434868data %>% filter(vacname2 == "Hep B vaccine for newborn") %>% nrow()## [1] 7737858Create a province vaccination coverage function based on the number of children having been recorded.
province_coverage_plot <- function(data_temp) {
province <- table(data_wide$province_reg, data_wide$year_dob)
province <- as.data.frame(province)
province <- province %>% rename(province_reg = Var1 , year_dob = Var2)
frequency_table <- table(data_temp$province_reg, data_temp$year_dob)
frequency_table <- as.data.frame(frequency_table)
frequency_table <- frequency_table %>% rename(province_reg = Var1 , year_dob = Var2)
province <- left_join(x = frequency_table, y = province, by = c("province_reg","year_dob" ))
province <- province %>% mutate(coverage = (Freq.x/Freq.y)*100)
province <- province %>% mutate(coverage = format(round(coverage, 1), nsmall = 1, scientific = FALSE))
province_wide <- province %>% select(province_reg, year_dob, coverage) %>% pivot_wider(id_cols = province_reg, names_from = year_dob, values_from = coverage)
gg <- ggplot(province, aes(x = factor(year_dob), y = as.numeric(coverage), group = factor(province_reg), color = factor(province_reg))) +
geom_line() +
labs(title = "Coverage by Year",
x = "Year",
y = "Coverage") +
facet_wrap(~province_reg) +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1),
axis.text.y = element_text(size = 5)) +
scale_y_continuous(breaks = seq(0, max(province$coverage), by = 20))
# Convert ggplot to an interactive plot using plotly
ggplotly(gg)
}The birth numbers for each region are based on number of children having been recorded.
region_coverage_plot <- function(data_temp) {
region <- table(data_wide$region, data_wide$year_dob)
region <- as.data.frame(region)
region <- region %>% rename(region = Var1, year_dob = Var2)
frequency_table <- table(data_temp$region, data_temp$year_dob)
frequency_table <- as.data.frame(frequency_table)
frequency_table <- frequency_table %>% rename(region = Var1, year_dob = Var2)
region <- left_join(x = frequency_table, y = region, by = c("region", "year_dob"))
region <- region %>% mutate(coverage = (Freq.x / Freq.y) * 100)
region <- region %>% mutate(coverage = format(round(coverage, 1), nsmall = 1, scientific = FALSE))
region_wide <- region %>% select(region, year_dob, coverage) %>% pivot_wider(id_cols = region, names_from = year_dob, values_from = coverage)
gg <- ggplot(region, aes(x = factor(year_dob), y = as.numeric(coverage), group = factor(region), color = factor(region))) +
geom_line() + # Use geom_line instead of geom_col for a line graph
labs(title = "Coverage by Year",
x = "Year",
y = "Coverage") +
facet_wrap(~region) +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1),
axis.text.y = element_text(size = 5)) +
scale_y_continuous(breaks = seq(0, max(region$coverage), by = 20))
# Convert ggplot to interactive plot using plotly
ggplotly(gg)
}data <- data %>%
arrange(pid, vacdate) %>%
group_by(pid) %>%
mutate(vac_order = row_number())
setDT(data)
data <- data %>% mutate(vac_order = paste("V", vac_order, sep = ""))data %>% group_by(pid) %>% mutate ()## # A tibble: 40,421,349 × 20
## # Groups: pid [12,477,883]
## pid vacname vacdate vactype vacplace vacplace_type province_reg
## <chr> <chr> <date> <chr> <chr> <chr> <chr>
## 1 1010101201400… ENGERI… 2021-11-13 TCMR Phòng t… DV Hà Nội
## 2 1010101201400… Hepava… 2017-07-24 TCCD Trung t… DV Hà Nội
## 3 1010101201400… Hepava… 2018-04-21 TCCD Viện đà… DV Hà Nội
## 4 1010101201400… Hepava… 2016-12-04 TCCD Trung t… DV Hà Nội
## 5 1010101201400… ENGERI… 2017-08-13 TCCD Trung t… DV Hà Nội
## 6 1010101201400… ENGERI… 2017-09-16 TCCD Trung t… DV Hà Nội
## 7 1010101201400… ENGERI… 2020-08-23 TCMR Trung t… DV Hà Nội
## 8 1010101201400… ENGERI… 2020-10-10 TCMR Trung t… DV Hà Nội
## 9 1010101201400… ENGERI… 2014-09-18 TCCD Bệnh vi… DV Hà Nội
## 10 1010101201400… Hepava… 2018-04-28 TCCD Trung t… DV Hà Nội
## # ℹ 40,421,339 more rows
## # ℹ 13 more variables: district_reg <chr>, commune_reg <chr>, sex <chr>,
## # dob <date>, ethnic <chr>, year_dob <int>, year_vac <int>, vac_age <drtn>,
## # vac_month <dbl>, region <fct>, vacgroup <chr>, vacname2 <chr>,
## # vac_order <chr>data_wide <- data %>%
subset(select =c(pid, vac_order, vac_month, province_reg, year_dob, region)) %>%
pivot_wider(names_from = vac_order, values_from = vac_month)
data_wide <- data_wide %>%
mutate(
HBVnewborn_delay = ifelse(V1 == -2 , 0,
ifelse(V1 == -1 , 1,
ifelse(V1 == 0, 2, 3))))
data_wide <- data_wide %>%
mutate(
HBV1_delay = ifelse(V1 > 8, 2,
ifelse(V1 >2 , 1,
ifelse(V1 == 2 , 0,
ifelse(V1 == 1, -1,
ifelse(V2 > 8, 2,
ifelse(V2 > 2, 1,
ifelse(V2 == 2, 0, -1))))))))
data_wide <- data_wide %>%
mutate(
HBV2_delay = ifelse(V2 > 9, 2,
ifelse(V2 > 3 , 1,
ifelse(V2 == 3 , 0,
ifelse(V1 >= 1, -1,
ifelse(V3 > 9, 2,
ifelse( V3 > 3, 1,
ifelse(V3 == 3, 0, -1))))))))
data_wide <- data_wide %>%
mutate(
HBV3_delay = ifelse(V3 > 10, 2,
ifelse(V3 > 4 , 1,
ifelse(V3 == 4 , 0,
ifelse(V2 >= 3, -1,
ifelse( V4 > 10, 2,
ifelse(V4 > 4, 1,
ifelse(V4 == 4, 0, -1))))))))
data_wide %>% filter(HBV1_delay == 0, year_dob == 2018) %>% nrow()## [1] 695231data_wide %>% filter(HBV1_delay == 1, year_dob == 2018) %>% nrow()## [1] 745797data_wide %>% filter( HBV1_delay == 0, year_dob == 2017) %>% nrow()## [1] 1011876data_wide %>% filter( HBV1_delay == 1, year_dob == 2017) %>% nrow()## [1] 608054data_wide %>% filter( HBVnewborn_delay == 0, year_dob == 2018, province_reg == "Bạc Liêu") %>% nrow()## [1] 13313data_temp <- data %>% filter(vac_month == -2)
interactive_plot <- province_coverage_plot(data_temp)
interactive_plotdata_temp <- data %>% filter(vac_month %in% c(-2,-1))
interactive_plot <- province_coverage_plot(data_temp)
interactive_plotdata_temp <- data %>% filter(vac_month == 2)
interactive_plot <- province_coverage_plot(data_temp)
interactive_plotdata_temp <- data %>% filter(vac_month == 3)
interactive_plot <- province_coverage_plot(data_temp)
interactive_plotdata_temp <- data %>% filter(vac_month == 4)
interactive_plot <- province_coverage_plot(data_temp)
interactive_plotdata_temp <- data_wide %>% filter(HBVnewborn_delay == 0)
interactive_plot <- region_coverage_plot(data_temp)
interactive_plotdata_temp <- data_wide %>% filter(HBVnewborn_delay %in% c(0,1,2))
interactive_plot <- region_coverage_plot(data_temp)
interactive_plotdata_temp <- data_wide %>% filter(HBV1_delay %in% c(-1,0,1,2) )
interactive_plot <- region_coverage_plot(data_temp)
interactive_plotdata_temp <- data_wide %>% filter(HBV2_delay %in% c(-1,0,1,2) )
interactive_plot <- region_coverage_plot(data_temp)
interactive_plotdata_temp <- data_wide %>% filter(HBV3_delay %in% c(-1,0,1,2))
interactive_plot <- region_coverage_plot(data_temp)
interactive_plotIn this report, we highlight a concerning trend regarding the up-to-date vaccination coverage of children born in 2019. According to a source, the coverage is low, primarily due to the impact of COVID-19. The current situation emphasizes the importance of addressing vaccination challenges, especially in the context of the ongoing pandemic.
data_temp <- data_wide %>% filter(HBVnewborn_delay == 0 )
interactive_plot <- province_coverage_plot(data_temp)
interactive_plotdata_temp <- data_wide %>% filter(HBVnewborn_delay %in% c(0,1,2) )
interactive_plot <- province_coverage_plot(data_temp)
interactive_plotdata_temp <- data_wide %>% filter(HBV1_delay %in% c(-1,0,1,2))
interactive_plot <- province_coverage_plot(data_temp)
interactive_plotdata_temp <- data_wide %>% filter(HBV2_delay %in% c(-1,0,1,2))
interactive_plot <- province_coverage_plot(data_temp)
interactive_plotdata_temp <- data_wide %>% filter(HBV3_delay %in% c(-1,0,1,2))
interactive_plot <- province_coverage_plot(data_temp)
interactive_plotlibrary(dplyr)
create_interactive_plot <- function(data_wide, selected_region) {
data_temp <- data_wide %>%
filter(HBVnewborn_delay == 0, region == selected_region)
interactive_plot <- province_coverage_plot(data_temp)
return(interactive_plot)
}
# Example usage
regions_of_interest <- c("RRD", "NMM", "NCC", "CHL", "SE", "MKD")
for (region in regions_of_interest) {
plot_name <- paste0("interactive_plot_", region)
assign(plot_name, create_interactive_plot(data_wide, region))
}interactive_plot_RRDinteractive_plot_NMMinteractive_plot_NCCinteractive_plot_CHLinteractive_plot_SEinteractive_plot_MKDlibrary(dplyr)
create_interactive_plot <- function(data_wide, selected_region) {
data_temp <- data_wide %>%
filter(HBVnewborn_delay %in% c(0,1,2), region == selected_region)
interactive_plot <- province_coverage_plot(data_temp)
return(interactive_plot)
}
# Example usage
regions_of_interest <- c("RRD", "NMM", "NCC", "CHL", "SE", "MKD")
for (region in regions_of_interest) {
plot_name <- paste0("interactive_plot_", region)
assign(plot_name, create_interactive_plot(data_wide, region))
}interactive_plot_RRDinteractive_plot_NMMinteractive_plot_NCCinteractive_plot_CHLinteractive_plot_SEinteractive_plot_MKDlibrary(dplyr)
create_interactive_plot <- function(data_wide, selected_region) {
data_temp <- data_wide %>%
filter(HBV1_delay %in% c(0,1,2), region == selected_region)
interactive_plot <- province_coverage_plot(data_temp)
return(interactive_plot)
}
# Example usage
regions_of_interest <- c("RRD", "NMM", "NCC", "CHL", "SE", "MKD")
for (region in regions_of_interest) {
plot_name <- paste0("interactive_plot_", region)
assign(plot_name, create_interactive_plot(data_wide, region))
}interactive_plot_RRDinteractive_plot_NMMinteractive_plot_NCCinteractive_plot_CHLinteractive_plot_SEinteractive_plot_MKDlibrary(dplyr)
create_interactive_plot <- function(data_wide, selected_region) {
data_temp <- data_wide %>%
filter(HBV2_delay %in% c(0,1,2), region == selected_region)
interactive_plot <- province_coverage_plot(data_temp)
return(interactive_plot)
}
# Example usage
regions_of_interest <- c("RRD", "NMM", "NCC", "CHL", "SE", "MKD")
for (region in regions_of_interest) {
plot_name <- paste0("interactive_plot_", region)
assign(plot_name, create_interactive_plot(data_wide, region))
}interactive_plot_RRDinteractive_plot_NMMinteractive_plot_NCCinteractive_plot_CHLinteractive_plot_SEinteractive_plot_MKDlibrary(dplyr)
create_interactive_plot <- function(data_wide, selected_region) {
data_temp <- data_wide %>%
filter(HBV3_delay %in% c(0,1,2), region == selected_region)
interactive_plot <- province_coverage_plot(data_temp)
return(interactive_plot)
}
# Example usage
regions_of_interest <- c("RRD", "NMM", "NCC", "CHL", "SE", "MKD")
for (region in regions_of_interest) {
plot_name <- paste0("interactive_plot_", region)
assign(plot_name, create_interactive_plot(data_wide, region))
}interactive_plot_RRDinteractive_plot_NMMinteractive_plot_NCCinteractive_plot_CHLinteractive_plot_SEinteractive_plot_MKD